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We describe a simple numerical simulation, suitable for an undergraduate project (or 
graduate problem set), of the Brownian motion of a particle in a Hooke-law potential well. 
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O ' Understanding this physical situation is a practical necessity in many experimental contexts, 

for instance in single molecule biophysics; and its simulation helps the student to appreciate 
the dynamical character of thermal equilibrium. We show that the simulation succeeds in 
capturing behavior seen in experimental data on tethered particle motion. 

Pq ■ PACS numbers: 05.20.-y, 05.40. Jv, 87.15. Aa 
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X> I. INTRODUCTION AND SUMMARY 

Introductory courses in statistical physics often place their greatest emphasis on average quan- 
ts- ■ tities measured in thermodynamic equilibrium. Indeed, the study of equilibrium gives us many 
powerful results without having to delve into too much technical detail. This simplicity stems in 
part from the fact that for such problems, we are not interested in time dependence (dynamics); 
accordingly dissipation constants like friction and viscosity do not enter the formulas. 

There are compelling reasons, however, to introduce students to genuinely dynamical aspects 
of thermal systems as early as possible — perhaps even before embarking on a detailed study of 
equilibrium^. One reason is that students can easily miss the crucial steps needed to go from a 
basic appreciation that "heat is motion" to understanding the Boltzmann distribution, and thus 
can end up with a blind spot in their understanding of the foundations of the subject. Although any 
kind of rigorous proof of this connection is out of place in a first course, nevertheless a demonstration 
of how it works in a sample calculation can cement the connection. 

A second reason to give extra attention to dynamical phenomena is the current surge in student 
interest in biological physics. Much current experimental work studies the molecular processes of 
life, or their analogs, at the single-molecule level, where simple mathematical descriptions do seem 
to capture the observed behavior. 

One familiar setting where simple models describe statistical dynamics well is the theory of 
the random walk, and its relation to diffusion. Ref. 1 showed an attempt to present classical 
statistical mechanics starting from the random walk, building on earlier classics such as Refi^. The 
link between thermal motion and the Boltzmann distribution emerges naturally in the analysis of 
sedimentation equilibrium: We require that in equilibrium, diffusive changes in a concentration 
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distribution must cancel changes caused by drift from a constant external force field (gravity). 

In this article we discuss a generalization of free Brownian motion that is important for inter- 
preting a large class of current experiments in single-molecule biophysics: The Brownian motion of 
a particle tethered to a Hookean spring. Specifically, we investigate the dynamics of fluctuations 
of such a particle in equilibrium. Although the theoretical analysis of this problem does appear in 
some undergraduate textbooks, it sometimes appears forbiddingly complex (e.g. Ref.^). We have 
found, however, that numerical simulation of the system is well within the range of an undergrad- 
uate project. Either as a project or a classroom demonstration, such a simulation brings insight 
into the emergence of equilibrium behavior from independent random steps, and also can serve as 
an entry into the topic of equilibrium fluctuations. 

The appendix gives an implementation written by one of us, an undergraduate at Ursinus 
College, in matlab. 

II. EXPERIMENTAL BACKGROUND 

As motivation, we briefly mention two contexts in which Brownian motion in a harmonic (Hooke- 
law) trap has played a role in recent biological physics experiments. 

Optical trapping-' 4 is now an everyday tool for the manipulation of micrometer-scale objects 
(typically a polystyrene bead), and indirectly of nanometer-scale objects attached to them (typ- 
ically DNA, RNA, or a protein). In this method, a tightly focused laser spot creates a restoring 
force tending to push a bead toward a particular point in space. When the trapping beam has a 
Gaussian profile, the resulting force on the bead is often to a good approximation linear in bead 
displacement. Thus the bead executes Brownian motion in a harmonic potential well. In such a 
well the motions along the three principal axes of the well are independent. 

The bead's motion in one or two dimensions can be tracked to high precision, for example by 
using interferomotry, thus yielding a time series. The probability distribution of the observed bead 
locations then reflects a compromise between the restoring force, pushing the bead to the origin, 
and thermal motion, tending to randomize its location. The outcome of this compromise is a 
Gaussian distribution of positions, from which we can read the strength of the harmonic restoring 
force ("trap stiffness"). For practical reasons, however, it is often more accurate to obtain both 
trap stiffness and the bead's effective friction constant from the autocorrelation function of the 
bead position (see Sect. IIVI below) . For example, slow microscope drift can spoil the observed 
probability distribution function (see Refi^). 

Our second example concerns tethered particle motion^. In this technique, a bead is physically 
attached to a "tether" consisting of a long single strand of DNA. The other end of the tether is 
anchored to a microscope slide and the resulting bead motion is observed. Changes in the bead's 
motion then reflect conformational changes in the tether, for example the binding of proteins to 
the DNA or the formation of a long-lived looped state. Fig. ^ shows example data for a situation 
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FIG. 1: Left, sample experimental data for the i-component of the motion of a bead of radius -Rbcad ~ 240 nm, 
attached to a DNA tether of length Ltcthor ~ 3500 basepairs, or ~ 1200 nm. The experiment, and the protocols 
used to remove drift from the raw data, are described in detail in Refp&. For clarity only the first 200 s of data are 
shown; the full data run lasted 600 s. Right, logarithm of the autocorrelation function of x expressed in nm 2 (see 
text). Solid line, experimental data. Dots, simulation described in this paper, using parameters A c g — 72 nm and 
effective viscosity 2.4 times that of water in bulk, [fitpm] 
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FIG. 2: Left, histogram of measured bead x position for the experimental data shown in Fig. Q The solid curve 
shows a Gaussian distribution with the same normalization and variance. Right, similar histogram for our simulation. 

[frhisto] 

where such conformational changes are absent, that is, simple tethered particle motion. 

As in the optical trap case, one can discard the dynamical information in the time series by 
making a histogram of particle locations. Fig. \Bp, shows the frequencies of occurence of various 
values of x. Rather detailed agreement between theory and experiment has been obtained for these 
histograms, including the slight deviation from Gaussian distribution shown in the figure^. In 
the present note, however, we are interested in a less sophisticated treatment of a more general 
question: Can we understand at least some aspects of the dynamical information contained in data 
like those in Fig. 

To this end, Fig.^> shows the logarithm of the autocorrelation function, C(r) = (x(t + r)x(t)), 
where the brackets denote averaging over t. At r = this quantity is just the meansquare dis- 
placement, which would diverge for a free particle but instead has a finite value determined by the 
equipartition theorem. At large times the autocorrelation falls to zero, because two independent 
measurements of x are as likely to lie on opposite sides of the tethering point as they are to lie on 
the same side. In fact, the autocorrelation function should fall exponentially with r, as we see it 
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does in Fig. [IJ>. 3 

III. SIMULATION SETUP 

We wish to simulate the motion of a bead of radius -Rbead) attached to a tether of length -Ltot) 
and compare our results to experimental data. For this we will need to know a specific property 
of DNA in typical solution: Its "persistence length" A is A ~ 45 nniji 

We will suppose that the external forces acting on the bead are a hard-wall repulsion from the 
microscope slide, a tension force from the tether, and random collisions with surrounding water 
molecules (see Sect. for further discussion). The tension force produces an effective potential 
well that keeps the bead close to its attachment point. In fact, at low relative extension, the 
tension exerted by a semiflexible polymer like DNA is approximately given by a Hooke-type lawi 
/ = — kx. The effective spring constant is k = 3/c B T/(2AL tot ), where k B T « 4.1 • 1(T 21 J is the 
thermal energy at room temperature and Ltot is the contour length of the polymer. Thus, again 
the motions in each of the x, y, and z directions are independent. Because the microscopy observes 
only x and y, we can reduce the problem to a two-dimensional one, and hence forget about the 
hard-wall force, which acts only along z. In fact, we will reduce it still further, by examining only 
the x coordinate of the bead position. 

There is a subtlety in that we do not directly observe the endpoint of the polymer in an 
experiment. Rather, we observe the image of the bead; the image analysis software reports the 
location of the bead center, a distance -Rbead from the attachment point. Thus the time series in 
Fig. reflects the motion of the endpoint of a composite object, a semiflexible polymer attached 
by a flexible link to a final, stiff segment of length -Rbead- To deal simply with this complication, we 
note that a semiflexible polymer can also be approximately regarded, for the purposes of finding 
its force-extension relation, as a chain of stiff segments of length 2A. In our case 2A ~ 100 nm is 
not much larger than -Rbead) so we will approximate the system as a single Hookean spring, with 
effective length L to t = ^tether + -Rbead and an effective persistence length A e g slightly larger than 
A. In the data we present, L t ether ~ 3500 basepairs, or « 1200 nm. We will fit the data to obtain 

Suppose we observe the bead at time intervals of At. Without the tether, the bead would take 
independent random steps, each a displacement drawn from a Gaussian distribution with mean- 
square step length 2D {At), where D is the bead's diffusion constant. If the bead were subjected to a 
constant force / (for example gravity), we could get its net motion by superimposing an additional 
deterministic drift on the random steps: Adrift^ = f/C- The friction constant £ is related to D 
by the Einstein and Stokes relations: £ = k^T / D = 67rr/i?bead; where rj is the viscosity of water. 
For the tethered case, at each step we instead use a position-dependent force — kx, where x is the 
current displacement. For small enough At (perhaps smaller than the actual video frame rate), x 
will be roughly constant during the step, justifying this substitution. 
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Here again we find a subtlety: The presence of the nearby wall creates additional hydrodynamic 
drag on the bead^. Moreover, the tether itself incurs significant hydrodynamic drag impeding the 
system's motion. Again for simplicity we acknowledge these complications by fitting to obtain an 
effective viscosity r] e g, which we expect to be larger than the value 10 -3 Pas appropriate to water 
in bulk. 

IV. SIMULATION RESULTS [s=sr] 

To summarize, the simulation implements a Markov process. Each step is the sum of a random, 
diffusive component with meansquare 2D(At), and a drift component —DKxjk^T. The constants 
D = k^T / (6Tn] c f{ i?bcad) an d n = Sk-gT / (2A c r L tot ) contain two unknown fit parameters, the ef- 
fective persistence length A e s and viscosity rj e s. The output of the simulation is the probability 
distribution of positions, and the autocorrelation function, which may be compared to experimental 
data. 

The simulation is deemed successful to the extent that the two fit parameters take values 
reasonably close to the expected values, differing in the expected directions, and the full functional 
forms of the outputs agree with experimental data. Figures Fig. © an d Fig. [21 show that indeed 
the simple model works well. Our simulation took At = 0.625 ms, for a total of about a million 
steps, which were then sampled every 40 ms to resemble the experimental data. 

V. DISCUSSION[s:e>] 

Our mathematical model made some naive simplifications. Two which have been mentioned 
involve the role of the bead radius, and certain sources of drag on the bead. In addition, there 
is a time scale for rearrangements of the DNA needed to change its extension, and for rotatory 
diffusion of the bead, which changes the location of the attachment point relative to its center. All 
of these effects have been assumed to be lumped into effective values of the fit parameters. 

Despite these simplifications, however, we did obtain two key qualitative aspects of the exper- 
imental data as outputs from the model, not set by hand. First, the autocorrelation function of 
equilibrium fluctuations has the expected exponential form. Moreover, as a consistency check, 
the experimental histogram of bead positions had roughly the Gaussian form we would expect for 
motion in a Hookean potential well. Both of these results emerge as statistical properties of a large 
number of simple steps, each involving only a diffusive step combined with a drift step based on 
the current bead location. 

The insights obtained from this exercise are different from those obtained from the analytical 
solution; for example, students can see the average behavior emerging from the random noise as 
the simulation size grows. In addition, the simulation approach opens the door to replacing the 
assumption of a harmonic potential by any other functional form. 
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APPENDIX A: MATLAB CODE 

function TetheredParticleAnalysis (Xdata,Ncorr,Nbins,deltat) 

"/This function simulates and analyzes the motion of a ID random walk 

"/confined in a harmonic potential well. 

% 

"/Written by Luke Sullivan, Ursinus College 

"/.Edited by John F. Beausang, University of Pennsylvania 

"/ 

"/Xdata = array of bead position (nm) 

"/nCorr = number of points in correlation function 

"/Nbins = number of histogram bins 

"/deltat = time step of data and simulation (sec) 

"/"/"/"/Experimental Data"/"/"/"/ 
Xdata = transpose (Xdata) ; 

n = length (Xdata) ; "/number of data points 

time = (1 :n) *deltat ; "/time series 

[FData,rData,histoData] = GaussHistoX(Xdata,Nbins) ; 
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logACData = LogAutoCorr (Xdata,Ncorr ,deltat) ; 
"/°/.°/.°/.S imul at ed dat a"/"/"/"/ 

Xsim = RandWalkSim (length (Xdata) ,deltat) ; "/generate simulated data 
[FSim,rSim,histoSim] = GaussHistoX (Xsim, Nb ins) ; 
logACSim = LogAutoCorr (Xsim,Ncorr ,deltat) ; 

'/.'/.'/.'/.output /. /;/;/. 

subplot (2,3,1) ;plot (time , Xdata, 'b' ) "/Experimental data 

title ('Experimental data'); 

xlabel ('Time (sec)'); 

ylabel ('Bead Position (nm) ' ) ; 

subplot (2, 3, 2) ;plot(rData,FData, 'bx' ,histoData(l, : ) ,histoData(2 , : ) , 'bo ') "/Gaussian curve 
title ('Gaussian fit (x) and Histogram (o) of exp data'); 
xlabel ('position (nm) ' ) ; 
ylabel ('Probability (1/nm)'); 

subplot (2, 3, 4) ; plot (time, Xsim, 'r')°/Simulated data 

title ('Simulated data'); 

xlabel ('Time (sec)'); 

ylabel ('Bead Position (nm) ' ) ; 

subplot (2, 3, 5) ;plot(rSim,FSim, 'rx' ,histoSim(l, :) ,histoSim(2, :) , 'ro') "/Gaussian curve 
title ('Gaussian fit (x) and Histogram (o)of Sim data'); 
xlabel ('position (nm) ' ) ; 
ylabel ('Probability (1/nm)'); 

subplot (2, 3, 3) ;plot(logACData(l, :) , logACData(2 , :) , 'b-' ,logACSim(l, :) ,logACSim(2, :) , 'r-')"/plot of correL 
title ('Compare Autocorrelation'); 
xlabel ('time difference (sec)') 
ylabel ( 'Log [Autocorrelation] , (nm2)') 

"/"/."/"/Subroutines"/"/"/;/"/ 

function [F,r ,Xhisto] =GaussHistoX (Xdata, Nbins) 

"/This function histograms the data and fits a Gaussian distribution 
Xmax = max (abs (Xdata) ) ; "/maximum position 

n = length (Xdata) ; "/number of data points 

binWidth = Xmax/Nbins; "/histogram bin width 

stdevX = std (Xdata) ; "/standard deviation of data 

F=zeros(l,n) ;r=zeros(l,n) ;Xhisto=zeros(2,Nbins+l) ; "/initialize 
Xhisto(l,:)= ( (1 :Nbins+l) - . 5) *binWidth; "/midpoint of histogram bins 

for i=l:n 

r(i)=abs(Xdata(i)) ; 

F(i)=2/sqrt(2*pi*stdevX~2)*exp(-(Xdata(i)~2)/(2*stdevX~2)) ;"/l sided gaussian curve 
which=l+f loor (abs (Xdata(i)) /binWidth) ; "/which bin data falls into 
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Xhisto(2,which)=Xhisto(2,which)+l; "/increment bin 

end 

Xhisto(2, : )=Xhisto(2, : )/n/binWidth; "/.convert counts to probability (1/nm) 
function logac = LogAutoCorr (Xdata,Ncorr ,deltat) 

"/This function determines the autocorrelation of the data for Ncorr points 
n = length (Xdata) ; 

logac = zeros (2, Ncorr) ; 

logac(l,:)= (0:Ncorr-l)*deltat; "/time steps 

for s = 1: Ncorr 

temp = zeros (l,n-s+l); 

for i=l : (n-s+1) 

temp(i)=Xdata(i)*Xdata(i+s-l) ; 

end 

logac (2 , s) =loglO (sum (temp) / (n-s+1) ) ; 

end 

function Xsim=RandWalkSim(n,deltat) 

"/This function simulates a ID random walk in a harmonic potential 
°/°/°/°/phy sical parameters "/ 7,7,7, 



L = 3477*. 34; "/tether length (nm) 

Xi =72; "/tether persistence length (nm) 

Rb = 240; "/bead radius (nm) 

kbT = 4.1*10~(-21) ; "/thermal energy (J) 

eta = 2.4*10" (-30) ; "/viscosity of H2D (J*s/nm~3) 

D = kbT/ (6*pi*eta*Rb) ; "/Stokes diffusion constant (nm2/s) 

kappa = 3/2*kbT/Xi/L; "/spring constant (J/nm2) 

mu = D*kappa/kbT*deltat;°/bias of step in harmonic potential 

ldiff = sqrt(2*D*deltat) ; "/diffusion length (nm) 



Xsim=zeros (1 ,n) ; 
for i=2:n 

deltaX=randn(l) *ldif f-Xsim(i-l) *mu; "/step size for each element called 
Xsim(i)=Xsim(i-l)+deltaX;7 new element value 

end 



